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Sh.S.Nickolajshvili, V,P,Polivanskii» 

At present Monte Carlo method is a basic one for the 
solution of the radiation transport problems in inhomogene- 
ous (finite and Infinite) media under conditions of energy 
dependence of the cross sections and in the process of ani- 
sotropic scattering* The main advantage of Monte Carlo me- 
thod in comparison with the generally used nvimerical metho- 
ds consists in suitability to multidimensional problems and 
bristly reveals in the transport theory problems. 

In fact, the solution of a general kinetic equation 
for the fluz<^ (Z'jJljE ) depended upon six variables appears 
as insurmountable task. In many cases when it is required the 
calculation of the functionals of a certain kind the problem 
of multidimensionality is not facilitated by the theoiy of 
adjoint functions since the dimension of the problem is 
still too high. 

At the same time according to the character of the in- 
formation necessary to receive Monte Carlo method may be 
easily modified for the purpose of automatic elimination of 
a nxunber of variables. Two variants of the method are deve- 
loped in a given report. One of them deals with the problem 

T ^ mm, 

of calculation of the flui<p ) for a certain rest- 

ricted set of coordinates, directions and energies, and the 
other is applicable to tbe problem of detemninatlon of the 
functionals of a certain kind. 

In many oases where radiation fields considered posn^a ..a 
tbe high degree of symmetry it is reasonable to use the dix- 
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ference methods. Thus, for example, the ntunerical solution 
of the problem for the axi^lar and spectral distributions 
of radiation at different distances from a point isotropic 
source is more effective in comparison with the well-known 
moments method. The region of application of difference 
methods may ba significantly expanded if to use the theoiy 
of adjoint equations. The latter is valid in that cases 
where the functionals sought for depend on a lesser number 
of independent variables in comparison with the function 
describing space— energy and angular distributioiis of radiar- 
tion. 

The difference methods. 

Plane-parallel problems. The solution of many problems on 
gammar-ray penetration through the material Is reduced to 
determination of the function C Z /i ) satisfying the 
equation 1 • 

( 1 ) 

CoJo(. C2) 

xmder supplementary conditions 


(jY(aAu./i)= yUCC? (3) 

The functions ^ , «x , 2. and G- are supposed known. Por the 
approximate solution of equation CO 1®^ us divide inter- 
val[ 06 z^:Qj into fl parts, intervalC'^ ^ i D into 
parts and along axis plott the points 
*...«^o is a least wave lenght of primary gamma^rays. The 
points of division of the intervals 2 C/jJ and 

4 denote asZic, K« 0,1, ...^n and 0,1,.*., 

,ni, JaaO, Q , /•*(>" respectively. 

Approximate values of the function vYCt ,^4,*^ ) in the point 
denote bb JYi,^ . 

Similar notation is used for the other functions. 

Let us return to the equation (O ®^d assume/-^ »A*j/ and 
A t Then we have 
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f 

Explanation is not required for notation. The function is 
defined by the integral (2) at ° and ^ a/lp. 

Suppose the right side of equation (4) is known, then 
the problem is reduced to the integration of ordinary first 
order differential equation. Integrating this equation by 
^ within and using the quadrature trapezoid 

rule receive approximately 

hAsrLrK )* 


z V -VKM -WK JUK /• 

^ ^k+r , Kb 0,1,,,,,,, n-1. Additional conditions 
following from the boundaiy ones (3} are necessary to add to 
the latter system 

( 6 ) 

After this, the system (5) is solved in the direction of 
increasing K ^ , and in opposite direction, if 

Proceeding to the approximate calculation of the fun- 
ction T (Z ,yU ) at the right of equation (4) note that 
its expression (2) representing usual record of the scat- 
tering integral is not suitable for our purposes. The me- 
thod proposed is based on Che other representation of this 
Integral obtained by the substitution of integration vari- 
ables and has the following form 
2jr 


o -t 


where 


+ cc?io( 


(7) 


( 8 ) 


It is supposed that 2. is not depend on Z • This 
restriction is not substantial. 

The division of the interval C”'* ^ lj] is perfor- 

med in that manner that none of the points is 
not fall on the point Al a 0, 
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Let/W * andj « Jp obtain 

where /,i' is calculated by formula (8) at /U ^ jUv • 

For calculation of the inner integral it may be constructed 
the quadrature formula with the points defined from 

the conditions 

^ = P-fS-ij - • .,P- (10) 

Ihen, write down the integral (9) approximately in 
iTp = JTnApp JY^ ^ ^ 

yd-|‘^‘codoC (12) 

^low, suppose that the values of the function IT(Z ) 

at vA =‘A,^<.Ap are already found. Then, by substituting the 
relation (11) into the equation (4) and transf ering 
into left obtain the following expression; 

^IV-P p ^ 

^ JjV(z,f{p^,ylp-<v) do< (13) 

The right side of the latter equation contains known func- 
tions and the problem is solved by the method decribed abo- 
ve. It remains to find the following integrals. 

Jp<^ = jjV'(z,pi,^p-^J C?|c^ ( 14 ) 

at M (See formula 12), For this purpose interval 

[]-1£^j6 + lJ is divided into two parts (—1,0) and (0, + l) and 
on each let be app:^oxlmately 


aon lei; oe app:^oximateiy 

( >(<o. 

, Z-i and 2, are some constants ai 


Here, Z± and 2^ a:^^ some constants and the coefficients 

Q and ^ may be determined by different ways, for exam- 

O ® XT 

pie, by the least- square method. 

By substituting the expression (15) into (14) obtain the 
integrals which are not difficult to calculate. Thus, the 
determination of space and angular gamma-ray distributions 
for every new value of energy variable is reduced to solu- 
tion of the ordinary differential equation (13) with known 
ri^ht side defined at previous steps of estimation. There- 
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fore the necessity of application of the successive appro- 
ximation method used for the difference solution of multi- 
velocity kinetic equations (2,3) is eliminated. It is also 
necessary to note that the decribed method may be used for 
the solution of neutron transport problems in light mode- 
rators without siguiflcant modifications. 


Spherically- symmetric nroblems. In this case the problem 
is reduced to the solution of the equation (2) 

IJ -f _i„ y // “Ti J 

A'sT" fl6) 

with boundary condition 


= ( 17 ) 

In the equation (16) the scattering integral 5-" is determi- 
ned by the same formula as in the case of the plane-parallel 
problems (with the difference that z is replaced byX), 
Hence, the problem is reduced to the differential eqxxation 
(16) with the known right side which may be solved by diffe- 
rence methods described in reports (3»4), 


Monte Carlo Methods. 
Calculation of the flux ^ (Z ), 


It is known that the method of direct modelling of par- 
ticle motion in a substance is suitable for the comparatively 
narrow class of problems, A number of modifications of the 
method was developed for increasing its statistic efficiency. 
One of this modifications of the method applicable to 
the radiation transpoit problem is a method cf local calcu- 
lation of the fltix. In the primary formulation cf the method 
(5>6) it is made use of the fact that the scattered part of 


the f lux <^ ( Z, , JL , £ ) in a fixed point of space Z*) 

is equal to th^ average of raudom variable 

(18) 

where indicate the position of the -th particle 


scattering, direction and energy after scattering, respecti- 
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( is a scattering function, i.e,, 

calculated for the solid angle unity the probability of a 
particle with the energy before scattering E n-i and directi- 
on -^w-i to scatter in the direction JT-n (it possess the 
energy E,^ a'f (En-i In the formula ( 18 ) W'n is 

a "weight" of the h -th scattering, T C'Z> , , E„) is an op- 

— ^ jk 

tical distance between Zn and Z points for the particles 
with energy E’n ,n (a^n. x) is a characteristic fxmction of 
the interval AXnear. x. 

Every individual realization of the process of parti- 
cle wandering ("history") gives the sampling values of par- 
ticle coordinates in a phase range Q(x ) by which 

may be calculated the corresponding value of random variable 
^ • According to the theorem proposed by A*J.Hinchin (?) 
the convergence on probability of the average from K "His- 
tories" towards the accurate value of the flux at N— 
takes place* Due to the fact that the higher moments of the 
random variable if and, in particular, the standard deviation 
do not exist this convergence is highly slow* The other 
deficiency of this method is in the necessity of survey of 
a very great number of histories for the obtaining of a 
reliable information about angular and energy distribution 
of the flux* 

M.Kalos (8) have proposed the modification of the me- 
thod mentioned above having the finite standard deviation 
of the total flux 

=SJ<p(Z*Ji.E)clJldE ^ 

This method is based on the inclusion of the multiplier 148,, 

causing the divergence of the standard deviation into 
the transition probability density of Markov’s chain of 
particle positions in a phase space and, therefore, changes 
the sampling character* 

It may be proposed the modification of Monte Carlo me- 
thod giving the standard deviation for the total flux and 
not changing the sampling character. The other valuable 
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feature of this method is a possibility of obtaining of 
energy distribution ,/l , £ ) for the separated fixed 

(iii*6c*t ion-iu • 

In the meth-od described above the estimation of the 
flxix is based on its representation in terms of the sum of 
non— scattered radiation contributions the sources of which 
are the separated points of particle scattering. It is easy 
to show that the similar representation takes place if to 
summarize the contributions of a single- scattered radiation 


so that 




■ 2 ,, (^n,En J X 


^ Eh-i 


JLv\-r 

an 




Jl - 


(19) 


where 


JL ■ tn , 'L, = " t.«-JL in 


In the formula (19)<J>' £ ) is a contribution of 
a single— scattered radiation^ M.E, ±s a symbol of mathema- 
tical expectation. 

It may be shown that the estimation, following from 
the formula (19) has the finite deviation for the total 
flux. The standard deviation of the flux does not exist, 
hovrever, the difficulty connected with a negligible phase 
volume in a space of directions is el i m i nated. 

As an example, consider the case of the infinite homo- 
geneous medium with isotropic scattering. For the Isotropic 
source the fomnila (36) leads to the followiiig result; 

where ***' St' 


UZ _ 7 - 


Approved For Release 2009/08/26 : CIA-RDP88-00904R000100110011-7 





Approved For Release 2009/08/26 : CIA-RDP88-00904R000100110011-7 


In the table 1 the results of ten sampling each inclu- 
ding 200 histories for the case "0,75, 'tZt obtain- 

ed by two methods mentioned above are presented. 

The other method of increasing of the efficiency of the 
local calculation method of the flux is based on that fact 
that the random quantities 

aBd 

possess the finite dispersion, if X and T are random vec- 
tors with non-degenerated distributions in D emd S, respec- 
tively, Here, and are arbitraiy square summable in D 
functions (on S, respectively), D is a region of a three- 
dimensional space of the variables and is a 

region of a surface in the same space. 

This circumstance allows to calculate the Fourier coef- 
ficients of the flux function sampling the orthonormal with 
weight in D functions as function. 

If the weight function P ( is successfully sampled 
the restoration of the flux function by its Fourier coeffi- 
cients is performed with a sufficiently high precision. Si- 
milar considerations are valid for the flux calculation at 
the surface S, 

The calculations are performed In the following manner. 
In the beginning of each”history" of a particle randomly InD 
it is sampled the position of the detector Z • The value 
^ ( Z*) is calculated. Further, the value ^ ( Z- ) may be 
additionally multiplied by^Ct*) and PCz*), _p 
vector uniformly distributed in D, then M.O,^ .( Z 
p( ) is a Fourier coefficient of the flux function, "fi (. "Z ") 
are orthonoimal in D with weight P( Z. ) functions (l»0,1,,,., 
n). 

However, this method of calculation of the Fourier coef- 
ficients is not the best one because H.O* J.( Z. )'^(L(Z. )P 
( may be a little quantity, while the dispersion 

C (■Z*)J may appear significant. In these 

eases the following generalisation of theorems 1 and 2 of 
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C 9 J is useful* Further, we use the notation of this report. 
Let (Q) he a square summable with weight P(Q)(P(Q) O) in 
D function (q6 d). 

Table 1, 


Flux 

Ho, samolin^ 

- 

C36y 

1 

0.0678 

0.0656 

2 

0.0689 

0,0660 

3 

0.0624 

0.0634 

4 

0.0843 

0.0721 

5 

0.0596 

0.0637 

6 

0.0688 

0.0628 

7 

0,0842 

0.0689 

8 

0.0746 

0.0660 

9 

0.0896 

0.0761 

10 

0.0567 

0.0615 

Average 

Exact value 

0.0709 

0.0665 

0.0666 


Consider a generalized interpolating polynomial 
P((9j =C<,f, ((P)+Cifi (Cpj ('(p; , constants Ci (i"0,1 . .,n) 

are defined from the conditions 

^ P«?i) 

where (Vj are certain fixed points of a region D, The follo- 
wing theorem Is valid. 

Theorem 1 * If the points Q are randomly sampled in 
D with the density of a probability Ww j ••• y Qn )» 

tnen 

M-O.Ci (Q) fi, (Q]c{c?, l = o.i,.. . . n. 

Under the same assumptions it is valid. 

Theorem 2, A standard deviation of the random value Ci. 

is equal to ^ ^ 

sS P ^ M-O-Cb o(^. 

u i.»0 
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results of the oalculation of a numerical gamma- 
spectrum in the air the local calculation method are pre- 
sented in the Table 2, 

The similar values obtained by a linear algebraic in- 
terpolation by random points are presented in the Table 3. 

In a given case the weight function is not introduced. 

The four groups of data are presented in these tables. 
Each group is obtained as a result of 1,000 tests, 

A point isotropic monoenergetic source with Eoa1,25HeT 
is used as a gamma— ray source. The source-detector distance 
is 30 meters, 0-1,25Mev energy Interval is divided into 
eight uniform groups. If the quantity presented in the tab- 
les 2 and 3 multiplied by 10 then we obtain a number of 
gamma- rays /second/cm , which is belonged to the given energy 
group. 

Table 2. 


No, of sampling 
No, of grouT) 

1 

2 

3 

4 

1 

1,580 

0,925 

0.737 

0.808 

2 

1,039 

1,162 

1.123 

0.975 

3 

0,450 

0,296 

0.342 

0.496 

4 

0,826 

0,303 

0.104 

0,228 

5 

0,123 

0.105 

0.214 

0.149 

6 

0,138 

0.594 

0.376 

0.129 

7 

0.053 

0,145 

0.110 

0.112 

8 

0.073 

0,208 

0.164 

0.022 

Integrated flux 

A. 832 

3.322 

4.1o7 

2.917 




Table 


1 

1,006 

0.864 

0,678 

0.827 

2 

1.463 

1.359 

0,910 

1,322 

3 

0.335 

0,437 

0.753 

0.173 

4 

0.301 

0,127 

0,117 

0.292 

5 

0.461 

0.415 

0.289 

0.576 

6 

0,398 

0.095 

0.281 

0.356 

7 

0.192 

0.050 

0.694 

0.000 

8 

0,000 

0.417 

0,087 

0,460 

Integrated flux 



3.809 
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V. 


Calculation of the functionals from the flux 

^ il .X.B'j, 

1. Monte Carlo method is often used for calculation of 
the functionals J- from the f lux ^ (l ,iL,f ) in a space- 
energy region in the inhomogeneous problems of a neutron 
transport theory 

iEo iVc AE AV C 20 ) 

vhere (p ; %$j^jiE^') is a Green function of the cor- 

responding transport equation* As a rule the regions of de- 
termination of the source S and the function 

(z ,ii ,E ) are greatly differed, therefore, the increase 
of a statistic efficiency is achieved when a sampling of the 
original coordinates of history is performed in accordance 
with a adjoint function in relation to functional^ CiQ 3 • 
Accurate finding of a adjoint function is nor less complex 
problem than the original one, therefore for evaluation 
of the functional 4 it is naturally to try to construct 
such adjoint function. If this ftmction Q ('Ic-fl.j £ 0 ) is 
constructed, the fxmctional 4 written in a following 

manner [lj ^ ^ ^ 

X JcIa. 4 Ej I ^i)~l > C21 ) 

where^ E^) q 

O / -y r- 1 J dJi JdE S- O 

K V.Za,-/lo,to) ' Q iitiK, E« ) 

and sampling of original coordinates of the history may he 
performed in accordance with density S and initial weight^. 

When determining the parameters of the resolved neutron 
levels according to the results of the radiation capture 
measurements the necessity is appeared to calculate the fun- 
ctionals of a type (20) which represent in this case the 
average absorption rate in a planar sample of thiolmess H 
of the neutrons from a plane source with a given spectrum 
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S(E^) in the energy interval 4 entirely comprising the 
resonance considered (11,12,13). In this problem Monte Car- 
lo method is used for evaluation of the average absorption 
rate of neutrons suffering scattering which has the 
following form: „ 

s (ej c/f E.Ej 

where ^ , £■,£,) is a neutron flux after a single scat- 

tering. The application of a ad;)oint function for estimati- 
on of a functional (22) is of importance as a region of 
determination of the source ^ always cojgipiderably exceeds 
the region where a adjoint 

is differed from zero. " 

For a construction of the approximate function C? (foj 
take a following model of the neutron transport. Let us be- 
lieve that the neutrons diffuse-through a sample of an effe - 
ctive thickness W without a collision but at the boundary 
some fraction of neutrons is reflected with an energy change 
by the magnitude of an average logarithmic loss ^ i 

Then an approximate adjoint function is written down 
in the following manner: 


Pm P n" St, 


CPi) 


1 fr- I 

where , r- (f^) , -nz(E/t)\ 

p . j_ 6k '^t^j J 


Gt (Ek) 




The normalized adjoint function is always used in an 
importance sampling of original coordinates, therefore the 
deviations of parameter values H*" and P from H and I, 
respectively do not influence greatly on the efficiency of 
approximate function application. 

The main feature of this function - a maximum displace- 
ment relative to a resonance energy — is well interpreted 
by introduction of an average logarithmic loss J into the 
expression (23). In the case of realization of Monte Carlo 
method for estimation of the functionals of a type (22) it 
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I 


/a -/j,/-. 


is important that introduction of importance sampling by a 
adjoint function (23) does not prevent from the application 
of the similar trajectoiy method for the simultaneous corre- 
lated calculation of a wide range of sample thicknesses. 

This method of improvement of a statistic efficiency 
is highly useful in a problem of calculation of a resonance 
absorption in a stratified activation detectors [‘14J, With 
the help of one approximate adjoint function for a sample 
of intermediate thickness it is managed to carry on the es- 
timation of activation absorption at a single level in the 
inner detector layer simultaneously for a wide range of 
filter thicknesses. These calculations were carried out for 
some activated isotopes with the purpose of obtaining some 
optimal sample thicknesses for the neutron spectrum measu- 
rements within a reactor. 

2, The other way of statistic efficiency improvement of 
Monte Carlo method in the case of estimation of functionals 
(20) when its region of determination A e xAV is conside- 
rably le^ ttan the region of determination of the flux den- 
sity f ) is based on obtaining of some information 

about an importance function in the process of sampling and 
its application for choice of further fate of a neutron, 

Ttiis is a method of splitting# 

If in the process of random sampling before next collisi- 
on a neutron has the coordinates not belonging to 

the region of a functional determination, then the probabili- 
ty to contribute to this functional after collision is 
equal to corresponding part Oj of a single volume of hyper- 
cube. Sampling randomly and xmlformly the point f from 
this part of a hypercube we obtain at every collision a 
contribution to functional and one of the branches for a 
further neutron fate. The second branch ensuring the com- 
pensation of a bias of such sampling corresponds to the 
random point from the hypercube region additional toG-;f . 
Rampling between these two branches is performed in a ran— 
dom^ manner in accordance with its weight factors Gj: and 

- 13 - 
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1— Gjj, respectively. 

The application of the method of splitting is highly 
useful in the case of Monte Carlo calculation of the influ- 
ence function for monoenergetic isotropic source located 
at the block surface in the infinite cylindrical lattice, 
Tlais fxmction being analogous to the Plachek’s function in 
a homogeneous medium is interest in the resonance absorption 
problems in heterogeneous media. 

In this problem the following quantity was considered 
as a functional j- . 

representing the average density of collisions in a block 
in the energy interval 

As a rule the determination of boundaries of a multidi- 
mensional region and sampling of random point from this 
region is a very difficult task, therefore, it is desirable 
to arrange the splitting by one of the variables. 

In the case considered a region ensuring the contribu- 
tion to the functional j is defined by the following condi- 
tion: i?r. 

where ^^ 4 , is a azimuthal angle of the neutron velocity direc- 
tion after scattering; is a block radius. The splitting 
is performed by one variable, jU , a cosine of scattering 
angle in a system of inertia. There, one can obtain the fol- 
lowing ratios for weights of the branches* 

k 1^0 ~ M 0^ 1 1 ^ - G 1 ; 

where _ 

while sinj-, 

p( COiOh) COiJn - CoSrA ccS 

M— is a nuclear mass. 

Calculation of a heterogeneous Plachek's function showed a 
considerable improvement of a statistic efficiency of Mon- 
te Carlo method at the expence of application of the split- 
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ting teclmiques particular in a ease of very thin blocks. 
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